    clc; clear;
fig6b_test(0.05)

function fig6b_test(k_input)


    % 参数矩阵定义
    alpha = 2.4;
    A = [-alpha, 0, 0;
          0, 0, 0;
          0, 0, 0];
    B = [1, -4, -3.5;
         0, 1, 2;
        -1, -4, 1.5];

    % phi0扫描范围，每隔0.2步长
    phi0_list = -0.5:0.05:0.2;

    % RK4积分参数
    dt = 0.01;
    T_total = 400; % 总时间
    steps = round(T_total / dt);

    % 初始状态固定
    x_init_base = [1e-6; 0; 0; 0; 0; 0; 0; 0]; % x11,x12,x13,x21,x22,x23,q11,q21
    % phi单独设置

    % 存储结果
    MLE_values = zeros(size(phi0_list));

    parfor i = 1:length(phi0_list)
        phi0 = phi0_list(i);

        % 构造完整初始状态9维，加入phi0
        x_init = [x_init_base; phi0];

        % 定义ODE函数句柄，封装参数
        f = @(t,x) mCNN_coupled_ODE(x, A, B, k_input);

        % 计算最大李雅普诺夫指数
        MLE = compute_max_LE(f, x_init, dt, T_total);

        MLE_values(i) = MLE;
    end

    % 绘图
    figure;
    plot(phi0_list, MLE_values, '-o','LineWidth',1.5);
    xlabel('\phi_0');
    ylabel('最大李雅普诺夫指数');
    title(sprintf('k = %.3f 时的最大李雅普诺夫指数 vs \phi_0', k_input));
    grid on;

end

% mCNN耦合ODE系统
function dx = mCNN_coupled_ODE(x, A, B, k)
    % x: 9维向量
    % x(1:3) = x11,x12,x13
    % x(4:6) = x21,x22,x23
    % x(7) = q11
    % x(8) = q21
    % x(9) = phi
    
    y1 = exp(-x(7)) .* x(1:3);
    y2 = exp(-x(8)) .* x(4:6);

    W = 1 - 2*exp(x(9));

    K_mat = diag([0, 0, k]);

    dx = zeros(9,1);
    dx(1:3) = -x(1:3) + A*y1 + B*x(1:3) + K_mat*W*(x(1:3) - x(4:6));
    dx(4:6) = -x(4:6) + A*y2 + B*x(4:6) - K_mat*W*(x(1:3) - x(4:6));
    dx(7) = x(1);
    dx(8) = x(4);
    dx(9) = x(3) - x(6);
end

% 计算最大李雅普诺夫指数的示例函数
function MLE = compute_max_LE(f, x0, dt, T_total)
    % 简易最大李雅普诺夫指数估计
    delta0 = 1e-7;
    steps = round(T_total/dt);
    trans_steps = round(0.8*steps); % 丢弃前80%暂态

    x1 = x0;
    x2 = x0 + delta0*randn(size(x0));

    sum_ln = 0;
    count = 0;

    for i = 1:steps
        % RK4积分单步
        k1_1 = f(0,x1);
        k1_2 = f(0,x2);

        k2_1 = f(0,x1+dt/2*k1_1);
        k2_2 = f(0,x2+dt/2*k1_2);

        k3_1 = f(0,x1+dt/2*k2_1);
        k3_2 = f(0,x2+dt/2*k2_2);

        k4_1 = f(0,x1+dt*k3_1);
        k4_2 = f(0,x2+dt*k3_2);

        x1 = x1 + dt/6*(k1_1 + 2*k2_1 + 2*k3_1 + k4_1);
        x2 = x2 + dt/6*(k1_2 + 2*k2_2 + 2*k3_2 + k4_2);

        dist = norm(x2 - x1);
        if dist == 0
            dist = 1e-16;
        end

        if i > trans_steps
            sum_ln = sum_ln + log(dist/delta0);
            count = count + 1;
        end

        % 重设x2微扰方向
        direction = (x2 - x1) / dist;
        x2 = x1 + delta0 * direction;
    end

    MLE = sum_ln / (count * dt);
end

